% Note: 

function [FlowVec,CoefMat,UnitTrProbMat]=Distribution_TrProbMat(PP,SS,Pir,DistPolicy)

%% Step 1: from Lagfp_Lagzo to f_zo (Re-evaluation Effect)
Lagfp           =   SS.DistApp.fz.State(:,1); 
LagID           =   SS.DistApp.fz.State(:,2);
N_State         =   size(Lagfp,1);

f               =   Lagfp.*(1)/(1+Pir);

%--------------------------------------------------------------------------
% Transition Matrix for each Dimension
%--------------------------------------------------------------------------
UnitTrMat       =   cell(2,1);

% to f
TempGrid_Fixed  =   SS.DistApp.TransFun.f.fromLogDiff(SS.DistApp.fz.UnitNode{1});
TempGrid        =   SS.DistApp.TransFun.f.fromLogDiff(f);
UnitTrMat{1}    =   DistApp_Hist_TrMatApprox(TempGrid_Fixed,TempGrid);
% to zoID
UnitTrMat{2}    =   SS.FunApp.ExoState.TrMat(:,LagID);

%--------------------------------------------------------------------------
% Combined Transition Matrix
%--------------------------------------------------------------------------
TrProbMat_1     =   DistApp_Hist_TrMatAssemble(UnitTrMat,SS.DistApp.fz.UnitNum(1:2));

%% Step 2: from f_z to fp_z (Consumption and Saving)
ID              =   SS.DistApp.fz.State(:,2);
N_State         =   size(ID,1);
%--------------------------------------------------------------------------
% Transition Matrix for each Dimension
%--------------------------------------------------------------------------
UnitTrMat       =   cell(2,1);

% to fp
TempGrid_Fixed  =   SS.DistApp.TransFun.f.fromLogDiff(SS.DistApp.fz.UnitNode{1});
TempGrid        =   SS.DistApp.TransFun.f.fromLogDiff(DistPolicy.V.fp);
UnitTrMat{1}    =   DistApp_Hist_TrMatApprox(TempGrid_Fixed,TempGrid);
% to zoID
UnitTrMat{2}    =   sparse(ID,(1:1:N_State)',ones(N_State,1),SS.DistApp.fz.UnitNum(2),N_State);
%--------------------------------------------------------------------------
% Combined Transition Matrix
%--------------------------------------------------------------------------
TrProbMat_2     =   DistApp_Hist_TrMatAssemble(UnitTrMat,SS.DistApp.fz.UnitNum(1:2));

%% Step 3: from Lagfp_LagID to fp_ID
CoefMat         =   (1-PP.XI)*TrProbMat_2*TrProbMat_1;
FlowVec         =   PP.XI*SS.DistApp.Entrant.QW;

UnitTrProbMat   =   struct('ReEval',TrProbMat_1,'ConSaving',TrProbMat_2);